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Abstract 

We analyze the influence of reflective boundary conditions on the 
statistics of Poisson-Kac diffusion processes, and specifically how they 
modify the Poissonian switching-time statistics. After addressing sim¬ 
ple cases such as diffusion in a channel, and the switching statistics in 
the presence of a polarization potential, we thoroughly study Poisson- 
Kac diffusion in fractal domains. Diffusion in fractal spaces highlights 
neatly how the modification in the switching-time statistics associated 
with reflections against a complex and fractal boundary induces new 
emergent features of Poisson-Kac diffusion leading to a transition from 
a regular behavior at shorter timescales to emerging anomalous diffu¬ 
sion properties controlled by walk dimensionality of the fractal set. 

Keywords: Diffusion, Poisson-Kac process, Hyperbolic stochastic models, 
Transport on Fractals, Emerging statistical properties. 

1 Introduction 

The first, and certainly most intuitive microscopic description of diffusion 
processes is certainly via particle Brownian motion, dating to a series of 
papers by A. Einstein [1J and J. Perrin [2]. Mathematically, the connection 
between this archetype of irreversibility (diffusion), and random motion at 
particle level can be formulated in a rigorous way, starting from the defini¬ 
tion of stochastic integration (Ito, Stratonovich, etc.) as a particular form of 
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Stieltjes integration, introducing the concept of stochastic differential equa¬ 
tions (Langevin equations), and deriving their statistical properties, thus 
finally obtaining a parabolic equation for the probability density function 
(the forward Fokker-Planck equation) [3j 0> 5]- The mathematical tool es¬ 
tablishing the connection between stochastic dynamics at microscales, and 
macroscopic diffusion processes is the concept of Wiener processes (also re¬ 
ferred to as mathematical Brownian motion). A one-dimensional Wiener 
process W(t) is a continuous stochastic process possessing independent in¬ 
crements that are distributed in a Gaussian way (with zero mean and square 
variance equal to the time elapsed). 

Therefore, a diffusion equation of the form 

d t p(x,t) = Ddlp(x,t) (1) 

can be viewed as the Eulerian description of a microscopic Lagrangian mo¬ 
tion described by the stochastic kinematics 

dx[t) = V2 D dw(t) , (2) 

where dw{t ) are the infinitesimal increments of a one-dimensional Wiener 
process in the time interval (t, t + dt). 

Both eqs. m and Q possess fundamental physical limitations. The 
diffusion equation (|Tj) is characterized by a unbounded propagation velocity 
that manifestly violates the basic principles of the theory of special relativ¬ 
ity. This property is one-to-one with the fractal nature of the graph of a 
realization x{t ) vs t of the stochastic differential equation ([2]), which, as well 
known, possesses fractal dimension 3/2 [B]. 

In order to solve the above mentioned problem, stochastic models pos¬ 
sessing finite propagation velocity have been proposed, essentially based 
on the use of Poisson processes for generating the stochastic perturbations 
(see Section ed b ® ® mu na. In the simplest case of one-dimensional 
spatial problems, the particle moves with constant and bounded velocity, 
and the stochastic perturbation acts simply by changing the velocity direc¬ 
tion (velocity switching), and the statistics of the switching times follows 
an exponential distribution. The kinematic equations associated with this 
simple random motion are relativistically consistent and describe a family 
of stochastic processes that almost everywhere (i.e. apart from the time 
instants when the switching of the velocity direction occurs) differentiable 
functions of time. 

The basic article addressing this simple stochastic model is due to M. Kac 
published in 1974 [8], in connection with a stochastic interpretation of the 
telegraphers’ equation. For the sake of historical correctness, the first paper 
dealing with this class of models is by S. Goldstein J7] that considered similar 
processes referred to as persistent random walk. Subsequently, Gaveau et al. 
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applied this model to derive, via analytic continuation, the one-dimensional 
Dirac equation of relativistic quantum mechanics [12j . 

A significant and fertile Literature originates from these observations, 
using Poisson processes as a model of bounded noise, or as a prototype of 
colored noise (due to the exponential decay of the correlation function), see 
Dans na cm na m, 19] just to quote some of the Literature covering 
different aspects of physical interest. 

Henceforth, we will refer to this class of processes as Poisson-Kac pro¬ 
cesses (see Section [2] for details). In the physical Literature, Poisson-Kac 
processes are often referred to as dichotomous noise. 

While at smaller timescales (smaller than the average time interval be¬ 
tween two consecutive velocity switching) the motion is smooth, at larger 
timescales Poisson-Kac processes retain emergent Brownian properties and, 
under certain assumptions on the value of the model parameters, they con¬ 
verge to the ordinary Brownian motion case (Kac limit) [8]. 

In a previous article [20], we have studied the limitations imposed by the 
wave-like nature of the Poisson-Kac processes on the representation of phys¬ 
ically plausible boundary conditions to be applied to the (Eulerian) balance 
equations for the probability density function (actually, for the probability 
partial waves, see further Section 0. 

In this article, we analyze the impact of reflective boundary conditions, 
i.e. of impearmability conditions in a closed system, on the statistics of the 
switching times of Poisson-Kac diffusion processes. 

Indeed, due to the wave-like nature of these processes, any form of reflec¬ 
tion modifies the exponential switching-time statistics. Reflection conditions 
act as an active modulation of the Poissonian statistics. Specifically, com¬ 
plex boundary conditions are able to modify the long-term properties of 
the Poisson-Kac diffusion processes determining anomalous diffusive scal¬ 
ings. This means that the long-term (emerging) fractal properties of the 
Poisson-Kac processes are altered by complex geometry of the boundary. 

The article is organized as follows. Section [2] provides a brief descrip¬ 
tion of Poisson-Kac processes functional to the further developments of the 
present article, by considering also two-dimensional models that are fur¬ 
ther analyzed in Section [5j Section [3] introduces the issue of the influence 
of reflection boundary conditions on the Poissonian statistics in a simple 
case, namely that of particles in a straight two-dimensional channel (a form 
of Knudsen diffusion problem [21]). Section [4] addresses another relatively 
simple, albeit interesting model, namely diffusion in the presence of a po¬ 
tential determining polarization effects at the walls [20] . Section [5] addresses 
the case of fractal boundary conditions by considering Poisson-Kac diffu¬ 
sion on fractal sets (treated in a continuum way, and not merely as a lattice 
diffusion problem) [221 [25] . This is the first application of Poisson-Kac pro¬ 
cesses in disordered systems. This model can be also viewed as a stochastic 
microscopic description of a linear viscoelastic transport problem (transport 
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problem with memory) in a fractal medium. However, the main focus in 
this article is on reflection conditions and how they modify the emerging 
long-term fractal properties of a Poisson-Kac diffusive trajectory. 

2 Statement of the problem 

Consider the one-dimensional problem of a particle moving in a deterministic 
potential U(x) under overdamped conditions. If p is the friction factor, the 
particle experiences a deterministic velocity field v(x) given by 

v(x) = -~d x U(x ), (3) 

V 

to which stochastic fluctuations are superimposed. Therefore, the equation 
of motion can be written as 

dx(t) = v(x(t)) dt + - dF stoca (t ), (4) 

V 

where F stoca ,(t ) is the stochastic force, and dF stoca (t ) its infinitesimal incre¬ 
ment. In the case of the classical overdamped Langevin equation, F stoC a,(t) = 
rj \J2 Dq w(t) where w(t) is a Wiener process [5]. In the case of the finite 
propagation velocity model proposed by Kac, F stoca (t ) is the integral over 
time 

^stoca(f) = br] [ (—1 ) x(r) dr, (5) 

Jo 

where y(f) is a Poisson process possessing the switching rate a > 0, i.e., 
E \x(t)/x{ 0) = 0] = at, and b > 0. 

The Poisson-Kac equation corresponding to eq. ([6]) in the presence of a 
deterministic velocity field v(x) is given by 

dx(t) = v(x(t)) dt + b (—l) x ^ dt . (6) 

Let us indicate with X(t) the stochastic process associated with the equation 
(Qj]l at time t. 

Since (—l) x ^) can attain solely two values ±1, the statistical descrip¬ 
tion of the process is based on the two partial probability density functions 
p ± (x,t), 

p ± (x, t)dx = Prob{X(f) € (x, x + dx), (— l) x ^ = ±1} , (7) 

satisfying the wave equations with recombination (dissipative wave equa¬ 
tions) 

d t p + (x,t) = —d x [{v(x) + b)p + (x,t)] - ap + (x,t) + ap~(x,t) 
d t p~(x,t) = -d x [(v(x) - b)p~(x,t)] + ap + (x,t) - ap~(x,t). (8) 

4 



As these equations correspond to forward and backward propagating waves, 
the partial probabilities p±(x,t) are also referred to as partial probability 
waves. 

The probability density function p(x,t) for X(t) at time t is the sum of 
the two partial probability densities p ± (x,t) 

p(x, t ) =p + (x,t)+ p- (x, t) . (9) 

Let us briefly address the properties of the stochastic trajectories asso¬ 
ciated with eq. (0, and the inclusion of impermeability conditions in finite 
domains. 

First, consider eqs. 0 and 0 in the free space i.e. as the unconstrained 
propagation along the real line. As \'(t) is distributed in a Poissonian way, 
the statistics of the switching times r, i.e. thta of the time intervals between 
two consecutive switchings of (—is distributed in an exponential way, 
according to the probability density functions gp(r ) 

g P (r) = ae~ aT . (10) 

This is the bare Poissonian result. The exponential probability density func¬ 
tion for the switchnig times gives rise to the linear recombination terms 
±(— ap + + ap ~) entering the balance equations for the partial probabilities 
0, determining the exchange from p + to p~ and viceversa per unit time at 
a constant rate equal to a. 

Correspondingly, if a and b are bounded, the graph of a generic real¬ 
ization of eq. 0 as a function of time t is an almost everywhere smooth 
function in all the open intervals between two consecutive switchings. If 
t = t* is a time instant at which a switching event occurs, solely the deriva¬ 
tive of x(t) is discontinuous at t = t*, still keeping a bounded left and right 
derivative. 

Let us suppose that eq. 0 describes particle motion in a box, and that 
x = 0 and x = L > 0 are the positions of the box walls, so that particles 
moving in (0, L ) cannot cross them or perform any sort of tunneling at 
them. Suppose for simplicity that the deterministic contribution v(x) is 
absent. Impermeability at x = 0, L implies total reflections for the partial 
waves p ± (x,t). This means that at x = 0 the backward propagating waves 
p~(x,t ) is totally reflected, and this determines the boundary condition for 
p + (x,t), i.e., 

P + (0,t) =p~(0,t) . (11) 

In a similar way, the forward wave p + (x, t ) is totally reflected at x = L, and 
this induces the complementary boundary condition 

P~(L,t) =p + (L,t) . (12) 

In terms of stochastic trajectories, reflection at an impermeable wall cor¬ 
responds to an additional switching of (—l) x w, externally induced by the 
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boundaries, that superimposes to the Poissonian exponential statistics. If 
this occurs at x = 0, this means that the wall forces the transition of (—l)xM 
from —1 to 1, and similarly at x = L, the wall-imposed switching from 1 to 
— 1 occurs. 

Consequently, the presence of impermeable boundary conditions pro¬ 
vides an active contribution to the statistical properties of the trajectories of 
the stochastic process ([HD determining a modification of the effective switch¬ 
ing time statistics g(r), which is no longer equal to the bare distribution 
gp(r) expressed by eq. (HOD . The analysis of g(r) and of its consequences 
on the emerging properties of the solutions of eq. ([6]) are the main issues of 
this article. 

To conclude this brief introduction to Poisson-Kac processes, let us men¬ 
tion the connection between eq. © and the classical Langevin equation 
driven by Wiener fluctuations. For a, b —>• oo, keeping fixed the ratio 
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the Poisson-Kac stochastic equation (|6]) converges to the classical Langevin 
equation 

dx(t) = v(x(t)) dt + V2D dw(t ), (14) 

and the probability density function p(x, t) converges to the solution of the 
Fokker-Planck equation (for mathematical details see also [23]) 


d t p(x, t ) = -d x [v(x) p(x, f)] = D d 2 p{x , t) . (15) 


This property is referred to as the Kac limit. 


2.1 Two-dimensional Poisson-Kac processes 

Since in Section[5]we consider Poisson-Kac processes in complex two-dimensional 
structures (fractal sets), it is useful to address here the extension of eq. ([6]) 
in higher dimensions, which is an issue scarsely addressed in the physical 
Literature [251 [26] , but very fertile in theoretical probability theory unun]. 

Although there are very interesting alternative extensions of Poisson-Kac 
processes in higher dimensions, mainly addressed by Kolesnik and coworkers 
mm , we focus here on the most natural and straightforward generaliza¬ 
tion of eq. & i n two-dimensional domains and in the absence of external 
potentials or deterministic velocity contributions. 

The simplest generalization of eq. ([6]) in two-dimensional spatial co¬ 
ordinates in the case of pure stochastic motion considers two independent 
Poisson processes Xi(t), X 2 (t), each of which acts on distinct spatial coordi¬ 
nates, i.e., 

dx(t) = b{-l) xi ^dt 

dy(t) = b(—l) X2 ^dt (16) 
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possessing the same statistical properties, i.e. 


E[xi{t)/Xi(0) = 0] = E[x 2 (t)/x 2(0) = 0] = at . (17) 


Analogously to the one-dimensional case, the statistical properties of the 
purely diffusive Poisson-Kac model (1161) are fully described by the system of 
four partial probability density function t), here x = (x,y), 

p( ±,:t )(x, t) dx = Prob |x(i) € (x, x + dx), (-l) Xl{t) = ±1, (-1)»W = ±ij ; 

( 18) 

where dx indicates the two-dimensional area element, and X(t) = ( X ( t),Y(t)) 
is the two-dimensional stochastic process described by eq. m at time t. 

For the sake of notational simplicity, let us use the notation p\ (x, t) = 
p(+’ +) (x,f),p 2 (x,f) =p ( "’ +) (x,f),p 3 (x,t) =p (_ " ) (x,t),p 4 (x,t) =p(+-)( x ,t), 
so that the overall probability density function for X(t)) is the sum of the 
Ph s, 
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P(x,i) = 5^Ph(x,i). (19) 

h =1 

It can be easily recognized that each p/j(x, f) describes the statistics of a 
particle ensemble moving with constant velocity b/ 3^, where 


Ph 



^2 (h — 1)7T + 7T \ 


, sm 


2 (h — l)-7r + 7 r 


( 20 ) 


h = 1,... ,4. As it regards the recombination dynamics, consider for exam¬ 
ple pi = p(+’+), and use the notation (±, ±) to indicate ((—l) Xl ^ = ±1,(—1) X2 ^ = ±l). 
Considering the statistics of an ensemble of particles, the switching in the 
unit time from (+,+) to (—,+) and (+, —) occurs with rate a The reverse 
processes possess obviously the same rate a. Consequently, the balance equa¬ 
tion for pi(x, t) is the Eulerian continuity equation for a particle ensemble 
moving with constant velocity b(3± which accounts for the above discussed 
recombination dynamics amongs partial waves, 


<9 t pi(x,f) = — 6 V • - 2api(x, t) + a(p 2 (x,f) + p 4 (x,i)) . (21) 


where V indicates the nabla operator with respect to the spatial coordinates 
x. In a similar way, the balance equations for the other partial waves follow 


<9tp 2 (x,f) = —6V • (/3 2 p 2 (x, t)) - 2ap 2 (x, t) + a(p 4 (x,f) + p 3 (x,f)) 
<9tp 3 (x,f) = -6V • (/3 3 p 3 (x,f)) - 2ap 3 (x, t) + a(p 2 (x,f) + p 4 (x,t)) (22) 
<9 f p 4 (x,f) = —6 V • ((3 4 p 4 (x,t)) - 2ap 4 (x,t) + a(p 4 (x,f) +p 3 (x,f)) . 

Gathering eqs. d2U)-(I22D and summing over the index “/t” of the partial 
waves, the balance equation for the overall probability density function fol¬ 
lows 

d t p(x, f) = V ■ J (x, t) , (23) 
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where J(x, t ) is the probability flux originating from the stochastic bivariate 
Poissonian perturbation 
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J(:M) = b ^2f3 h p h (x,t). (24) 

h =1 

From eqs. m-m, it is straightforward to obtain the constitutive equation 
for the probability flux, multiplying each evolution equation for ph(x,t) by 
b 1 3 h and summing over the index “h ”, 

d t J(x,t) = -b 2 V • — 2a J(x, t). (25) 

Several observations follows from the system of equations (I23D - (I25|) : 

• the statistical description of the process corresponds to that of a linear 
viscoelastic material with memory, as the constitutive equation for 
J (x, t) depends on the first-order time derivative of the flux itself; 

• this representation is however irreducible in terms of p(x, t ) and J(x, t ), 
as the correct description of the process involves the full system of 
partial probability waves {_p^(x, t)}^ =1 ; 

• although there are some analogy between eq. (1251) and the Catta- 
neo transport equation with memory in higher dimensions 1231281, 
the functional form of the constitutive equation (1251) marks the funda¬ 
mental difference between the two models. For finite values of a and 
b, it is impossible to express the time-derivative of the probability flux 
J (x, t) solely as a function of J(x, t) itself and of the overall probability 
density function p(x, t): all the systems of {p^(x, t)}^ =1 should be con¬ 
sidered. If the first term at the r.h.s. of eq. (1251) is approximated by a 
term proportional to the gradient of p(x, t), as in the two-dimensional 
Cattaneo model [23 .28], the stochastic interpretation of the resulting 
constitutive equation is completely lost, and this is the fundamental 
reason while the higher-dimensional Cattaneo model does not satisfy 
the requirement of positivity |2fl| : 

• the classical diffusion equation follows from (1231) - ()25|) in the Kac limit. 

To prove the last observation, consider the limit for arbitrarily large 
a and b (a, b —>• oo), keeping fixed the ratio b 2 /a to a constant value. If 
the switching rate tends to infinity, the recombination process between par¬ 
tial waves becomes infinitely fast, and the structure of the partial waves 
“thermalizes” at any x, meaning that all the p^(x, t) collapse into a unique 
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function that is just p(x,i)/4. Enforcing this property into eq. (l25l) the 
constitutive equation for the probability flux becomes 



(26) 


where B is the diadic tensor 


B - 4 ^PhPh. ■ 


(27) 


h= 1 


Because of the structure of the normalized velocities (3 h defined by eq. (1201) , 
this tensor is isotropic and equal to the identity tensor, B = I. In the Kac 
limit, the second term at the right hand side of eq. ()26li vanishes and the 
constitutive equation reduces to the Fickian form, J(x, t) = —DVp(x, t), 
where the diffusivity is given by eq. m- 

Finally, let us discuss the role of impermeability conditions. Consider 
for simplicity the problem of Poisson-Kac diffusion inside a square, (.x , y) E 
(0, L) x (0, L). Because of the symmetries of the normalized velocities (3 hl 
each reflection at the box walls, placed at x = 0, L, y = 0, L corresponds 
to an externally forced switching either of (— l) Xl M (for the reflections at 
x = 0, L) or of (—l)X 2 (*) (f or thg reflections at y = 0 ,L). The boundary 
conditions for the partial waves at the walls can be obtained either by en¬ 
forcing that the normal component of the probability flux should vanish at 
the box walls or by expressing total reflection in terms of the partial waves 
as it result from the switching of one of the (—i)x/i(*), h = 1, 2. For example, 
consider the collisions with the boundary at x = 0. In this case (—l)* 1 ^) 
switches from the value (—ipK*-) = _i to the value (—l) Xl C+) = an d 
consequently the boundary conditions become 


pi(0, y, t ) = p 3 ( 0, y, t ), p 4 (0, y, t) = p 2 ( 0, y, t ), 


(28) 


meaning that the two backward waves (along the x coordinate) p 3 , and p 2 
are reflected back into p± and p 4 , respectively. Analogous expressions can 
be derived for the other wall collisions. 

3 The simplest case: Knudsen effect 

The simplest physical situation where the effects of impermeable walls con¬ 
trols the statistics of Poisson-Kac diffusion processes is represented by par¬ 
ticle motion in straight channels under creeping conditions, so that the de¬ 
terministic velocity fields acts only along the axial direction. Consider a 
two-dimensional model letting x be the axial, and y the transverse coordi¬ 
nates. 
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Particle motion is described (using nondimensional variables) by the sys¬ 
tem of stochastic equations 

dx(t) = v(y(t)) dt + b (—1) X1 ^ dt 

dy(t) = b(-l) X2{t) dt (29) 

where y £ (0,1), equipped with reflection boundary conditions at y = (1.1 
which represent channel walls. As the transverse motion is completely de¬ 
coupled from the axial one, we can consider the second equation (f29l) as a 
stand-alone problem. The Poisson process \2 (t) is characterized by a switch¬ 
ing rate a, and we assume b and a to be related by the Kac condition (1131) . i.e. 
b 2 = 2 Da, where D is the nondimensional diffusion coefficient. Therefore, 
instead of a and b we can use b and D or a and D to define completely the 
characteristic of the stochastic perturbation that controls transverse particle 
motion, ppp Figure |T] depicts the time series of the consecutive switching 
intervals r n vs the numeral order n. The most evident effect of impermeable 
boundaries is that the switching intervals r n are upper bounded, as the max¬ 
imum attainable value of the switching time is given by r max = L/b, where 
L is the transverse width of the channel, and b the characteristic velocity of 
the stochastic perturbation. In the present case, L = 1, and b = \/2 D a, so 
that 

1 , , 

T max — 7 = ■ (30) 

V2 U a 

For D = a = 1, r max = l/\/2, corresponding to the horizontal dashed line 
in figure [U 

This result follows from the elementary observation that the longest 
switching time is reached by a particle hitting consecutively the two walls 
without switching the velocity direction due to the Poissonian statistics. 
Consequently, the probability density function g{r ) for the switching times 
admits a compact support represented by the interval [0,r max ], while the 
bare Poissonian gp(r) does not. 

This effect can be referred to as the Knudsen effect on the Poissonian 
statistics induced by wall reflection, to mark the analogy with the phe¬ 
nomenon of Knudsen diffusion for dilute gases in narrow pores where the 
collision with the pore walls dominate the diffusive transport process. 

Figure [2] depicts the behavior of the switching time pdf g(r) obtained 
from the stochastic simulation of eq. (1^1) . The pdf g{r) can be approxi¬ 
mated by the superposition of a smooth exponential distribution defined in 
[0, r max ] possessing a decay exponent a e ^ a, different from the Poissonian 
switching rate a, and of an impulsive contribution centered at r max , i.e., 

g(r) = Ae~ aeT + B 5(t - T ma , x ), (31) 

where A and B are positive constants, such that A(l — e _ “ eTmax )/a e + l? = 1. 
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Figure 1: Time series of the switching intervals r n vs n for a Poisson-Kac 
diffusive trajectory at D = 1, a = 1. 

The effective rate depends on a and D, and its rescaled behavior, i.e., 
(a e — a)/a, is depicted in figure [3] as a function of a for different values of 
D. As can be observed, the relative deviation from the Poissonian rate a is 
more pronounced as D and a increase. 

It is possible to obtain a quantification of the distortion effects induced 
by wall reflections using simple arguments. Consider the fraction </> wa u of 
the number of switching events associated with wall collisions with respect 
to the overall number of switchings. A lower bound for </> wa u can be ob¬ 
tained by considering that this fraction should necessarily be greater than 
the Poissonian probability of having switching times greater than r max , i.e., 

<?Wi > / gp{r ) dr = e _aTmax = e ~^ a,2D = 0 o (a, D). 

J Tmax 

Figure0]depicts the behavior of </> wa n and of its lower-bound estimate D) 
as a function of a for two values of D. The estimate 4>o(a,D) approaches 
better the value </> wa n at lower values of a and higher values of D. 



4 Polarization model 

As a second example, consider the one-dimensional model in the presence 
of a deterministic biasing field, 

dx(t) = v(x(t)) dt + b (—l) x ^ dt , (33) 
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Figure 2 : Switching time probability density function g(r) vs r for the one¬ 
dimensional Kac diffusion process on the unit interval for a = 1 at different 
values of the diffusivity D. Line (a): D = 1, line (b): D = 0.1, line (c) 
D = 0.03. Line (d) corresponds to the pure Poisson statistics gp(r) = ae~ aT 
in the absence of boundary condition effects. 


defined for x G [0,1], equipped with impermeability conditions at x = 0,1. 
The velocity vield v{x) is given by 

v(x) = -v m sin , (34) 

and v m = 7 b, with 7 € (0,1), and a and b are such that D = 1. This case has 
been analyzed elsewhere [ 20 ], as a model for describing polarization effects 
within the Poisson-Kac paradigm, arising as a consequence of a nonvanishing 
value of the deterministic velocity field at one of the boundary points (x = 1 ), 
determining the formation of a polarization layer in the neighbourhood of 
x = 1 . 

In this article, we use this model to investigate better the conditions to be 
imposed on the stochastic dynamics (1331) in the presence of a nonvanishing 
deterministic bias at the boundary. 

First, consider the partial waves associated with eq. (1351) . i.e. its Eulerian 
description. The overall probability density flux J tot (x,f) is the sum of the 
convective (deterministic) and stochastic contribution, 

Jtot(x,t) = v(x) [p + (x,t) +p~(x,t)] + b [p + (x,t)-p~(x,t)] . (35) 
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Figure 3: Rescaled effective switching rate ( a e —a)/a vs diffusivity D at three 
different values of a. Line (a) and (□): a = 1, line (b) and (o): a = 0.1, line 
(c) and (•): a = 0.01. 

Since u(0) = 0, the zero-flux condition at x = 0 simply becomes 

p+(0,t)=p-(0,t). (36) 

This is the total reflection condition. At x = 1, u(l) =7 b, and the zero-flux 
boundary condition becomes in terms of partial waves 

p~(l,t) = \±lp + (l,t). (37) 

1-7 

Eq. (1371) is valid for 7 < 1. The extension to 7 > 1 is discussed in [20], and 
therefore it is not repeated here. 

For the stochastic differential equation (1331) the conditions at the bound¬ 
ary are simply reflection conditions, corresponding to the mirror simmet- 
ric reflection of the particle position at x = 0,1, and to the switching of 
(_l)x(b —(_l)x0) 0 f the Poissonian perturbation. 

In order to check the validity of these conditions, we compare the sta¬ 
tionary solutions p±(x) of the associated partial wave model (jU) with v(x) 
given by eq. (|34[) . in the presence of the boundary conditions m-m, and 
the overall stationary probability density function p*(x) = pf(x) + p*{x), 
with the corresponding stationary distributions obtained by integrating the 
stochastic equation of motion (1331) with an integration step ht = 10 -4 for an 
ensemble of N = 10' particle in which the above mentioned simple reflection 
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Figure 4: <^> wa u vs the switching rate a for two values of D. Symbols (□, •) 
represent simulation results, dashed lines the lower bound eq. (13211 . Lines 
(a) and (b) and symbols (□): D = 1, lines (c) and (d) and symbols (•): 
D = 10. 

conditions have been used at the boundaries. This comparison is depicted 
in figure [5] and is fully satisfactory. 

In point of fact, a more stringent test on the influence of the bound¬ 
ary conditions (especially at x = 1 ) is based on the analysis of the ratio 
P*( x )/pt( x )i that, because of eqs. (1361) and (1371) should attains the values 
1 at x = 0, and (1 + 7 )/(l — 7 ) at x = 1. Figure [ 6 ] shows the estimate of 
this ratio obtained from the stochastic evolution of the particle ensemble 
compared with the corresponding quantity obtained from the partial wave 
model at different values of 7 . 

The dashed horizontal lines correspond to the limit values (l+ 7 )/(l — 7 ) 
at x = 1 for the values of 7 considered. Also in this more severe test, the 
stochastic model equipped with the simple reflecting boundary conditions 
provide excellent quantitative agreement with the Eulerian results. 

The switching time distributions g[r ) obtained from stochastic simula¬ 
tions is depicted in figure [7] for several values of 7 . Compared with the 
corresponding results obtained for the simple Knudsen model described in 
Section [3l there are qualitative analogies and differences. The analogies are: 
(i) the compactness of the support of g(r), and (ii) the occurrence of an 
impusive contribution centered at r max . The differences are: (i) a significant 
modification in the low-r region induced by the polarization effects of the 
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Figure 5: Stationary probability distribution p*(x) (lines (a) and (□)), and 
stationary partial waves pf(x) lines (b) and (o)), p~{x) (lines (a) and (•)) for 
the polarization model (1331) atD = l,a = 5, 7 = 0.9. Solid lines represent 
the results of the integration of the partial wave equations, symbols the 
results of stochastic simulations of eq. (1551) . 

velocity field near the boundaries, and (ii) the exponential “backbone” of 
this distribution, that apart for impulsive and nearly impulsive contribution 
possesses an exponential decay equal to a, i.e., equal to the bare Poissonian 
switching rate (compare the decay of the Poissonian statistics, line (e) with 
mean local decay of the other lines (a) to (d)). 

5 Kac diffusion on fractals 

As a final example, consider diffusion in complex geometries such as in two- 
dimensional connected sets possessing fractal (Hausdorff) dimension 1 < 
df < 2. It is known that Brownian diffusion in fractal media possesses 
anomalous properties [ 221122 ], typically expressed by the power-law scaling 
of the mean square displacement (r 2 (t)) of Brownian particles 

(r 2 (t)} ~ t 2/dw , (38) 

where d w > 2 is the walk dimension [d w = 2 in Euclidean media, giving rise 
to the classical Einsteinian linear relation of ( r 2 (t )) ~ t). 

Most of the existing numerical results on diffusion in fractal media refer 
to lattice simulations, where Brownian particles move in a random and un¬ 
correlated way from a site of the fractal set to one of its neighbouring sites 
belonging to the set itself. 
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Figure 6 : Ratio p~ (x)/p+(x) under stationary conditions at D = 1, a = 5. 
Solid lines represent the results of the integration of the partial wave equa¬ 
tions, symbols the results of stochastic simulations of eq. (1331) . Arrows 
indicate increasing values of 7 = 0.7, 0.8, 0.9. Dashed horizontal lines cor¬ 
responds to the limit value at x = 1 , (l+ 7 )/(l — 7 ). 


In order to give concrete examples, consider the two fractal structures 
depicted in figure [ 8 j the Sierpinski carpet, and a loopless deterministic 
fractal, henceforth referred to as the deterministic cross fractal. Their fractal 
dimensions are df = log 8 /log 3 ~ 1.893 (for the Sierpinski carpet), and 
df = log 5/log 3 ~ 1.465 for the deterministic cross fractal. Anomalous 
diffusion properties on the Sierpinski carpets have been addressed in [22> 30j. 

Figure [9] depicts the behavior of the mean square displacement of clas¬ 
sical Brownian particles obtained from discrete lattice simulations. The 
simulations refer to a prefractal approximation of the structures at the 7-th 
iteration of the obstruction process, that corresponds to lattices possessing 
3' = 2187 lattice size per Cartesian coordinate. In point of fact, figure [ 8 ] 
depicts just these prefractal approximations, as can be noticed by the labels 
of x and y coordinates. 

From the scaling of the mean square displacement (line a) we obtain the 
numerical value 2/d w = 0.94 for the Sierpinski carpet (no analytic results are 
available for the Sierpinski carpet, which is an infinitely ramified structure 
not easily amenable to exact real-space renormalization). In the case of the 
deterministic cross fractal, the scaling theory of loopless fractals, based on 
the concept of chemical distance, predicts a simple relation between d w and 
dj , namely [ 22 . 231 

d w = df + 1, (39) 
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Figure 7: Switching time probability density function g(r ) vs r for the 
polarization model (I33l) - (1M1) at a = 5. Line (a): 7 = 0.7, line (b): 7 = 0.8, 
line (c): 7 = 0.9, line (d) 7 = 0.95. Line (e) represents the pure Poissonian 
distribution gp(r) = ae~ ar . 


and the numerical results perfectly agree with the scaling (r 2 (f)) j2/(df+l) 

(line b). 

Next, consider Poisson-Kac diffusion in these sets. The mathematical 
setting of the problem has been developed in paragraph 12.11 both as it 
regards the stochastic equations of motion and their statistical (Eulerian) 
description in terms of partial waves. In the numerical simulations of eq. 
m we consider that the unit site of the prefractal approximation possesses 
a unit linear length. Consequently, the unit square is the building block 
of the fractal structure, that at iteration n = 7 extends over a length of 
L = 3 7 = 2187. Using eq. (I16I1 . and applying total reflection conditions 
whenever a particle hits the boundary (which implies also the switching of 
(—l) Xl d) or (— 1)* 2 W depending on the which boundary is involved), it is 
possible to perform an off-lattice (continuous) simulation of the Poisson-Kac 
diffusion process. Let D = 1 be the nondimensional diffusion coefficient. 
Figure [12] panels (a) to (c) depict a portion of an orbit of a Poisson-Kac 
particle diffusing inside the Sierpinski carpet at different values of a. An 
analogous picture for the deterministic cross fractal is shown in figure fill 
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Figure 8: Panel (a): Sierpinski carpet. Panel (b): deterministic cross fractal. 


At small values of a, figure flOl panel (b) and figure [11] (keeping D fixed 
and consequently, varying b according to the Kac condition eq. (11.311 1. the 
orbits are significantly different from that of classical Brownian motion, as 
the switching of the Poissonian processes is controlled by collisions with the 
boundary of the fractal structure. As a increases (see panel (b) and (c) 
in figure [Toll , particle orbits resemble more closely that of classical Brow¬ 
nian particles. Figure [12] depict the mean square displacement of Poisson- 
Kac particles in the two fractal structures considered at several values of 
the switching rate a. The mean square displacement ( r 2 (t )) possesses two 
asymptotic scalings: at short time scales the mean square displacement is 
ballistic, and this is to be expected from the wave nature of the stochastic 
process, as for t « min{l/6,1/a} particles move in straight lines. In the 
long-term limit, the mean square displacement possesses exactly the same 
anomalous behavior characterizing classical Brownian particles. 

Anomalous diffusion effects, i.e. the occurrence of a value of d w > 2, 
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Figure 9: Mean square displacement (r 2 (t)} vs t for classical Brownian lattice 
diffusion on the fractal structures considered. Symbols are the results on 
lattice random simulations: (o): Sierpinski carpet, (•): deterministic cross 
fractal. Lines (a) and (b) represent the scaling ( r 2 (t )) t 2 / dw , with 2 /d w = 
0.94 (line (a)), and 2/d w = 2/(df + 1), df = log 5/log3 (line (b)). Line (c) 
represents ( r 2 (t )) r\_/ ' 


are in the Poisson-Kac case the emerging feature of the collisions with the 
boundary of the fractal structure, and therefore it is to be expected that 
the switching time statistics becomes significantly modified. As the process 
is two-dimensional, we consider separately the statistics g a (r), a = 1,2 for 
the switchings of the two Poissonian processes. The hit of the boundaries 
orthogonal to the x-axis determines a switching of (—and the re¬ 
flection onto a boundary orthogonal to the y-axis determines a switching of 
(_l)X2(i). 

Figure [T3l and ITLl depict the behavior of g Q (r) a = 1,2 obtained from 
stochastic simulations. The two density functions g\ (r), g-iij) are equal and 
differ from the pure Poissonian statistics gp(r). Nevertheless, their behavior 
is well approximated by an exponential decay with an effective decay rate 
a e different from a. 

A coarse information of the role of the boundary collisions can be ob¬ 
tained by considering the ratio 4> C o\i of the overall number of collisions with 
the walls to the number of switchings in the velocity determined by the 
Poissonian statistics. The behavior of 0 co ii vs a is depicted in figure [J5] at 
D = 1. It can be observed that this ratio follows a power law behavior 
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Figure 10: Realizations of diffusive Kac orbits on the Sierpinski carpet at 
D = 1. Panel (a): a = 0.1, panel (b): a = 1, panel (c) a = 10. 
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Figure 11: Realization of a diffusive Kac orbit on the deterministic cross 
fractal at D = 1 and a = 0.1. 

^coiKo) ~ a ~^, where the exponent £ attains the value (from the best fit 
of simulation data): ( = 0.72 for the Sierpinski carpet, £ = 0.79 for the 
deterministic cross fractal. This power law scaling is nor related in a simple 
way to the characteristic geometric ( df ) and dynamic ( d w ) exponents of the 
medium. We leave a more careful theoretical analysis on the properties of 
this scaling law, and of its framing within the theory of transport in fractal 
media to future investigations. 

Let us analyze in more detail the trajectories of Poisson-Kac particles 
on fractals. By definition, the trajectories of the realizations of Poisson-Kac 
processes are with probability 1 almost eveywhere smooth curves. Disconti¬ 
nuities in the derivatives arise as a consequence of the “internal” Poissonian 
switching occurring in average at a characteristic time r c = 1/a, and as a 
result of the “external” reflections with the boundaries. Since the unit build¬ 
ing block of the fractal structure possesses unitary characteristic length, the 
average time between collisions with the boundary is t w = 1/b = l/\/2 Da. 

For time scales much smaller than the minimum between r c and t w , 
the trajectories are simply straight lines. Figure [16] depicts two examples 
of particle trajectory (actually its x-component) at two different values of 
a. The higher a is, the smaller is the characteristic timescale at which the 
realization of Poisson-Kac processes attain asymptotic fractal properties. 

The role of wall collisions that superimpose to the intrinsic Poissonian 
switchings is therefore to induce in the long-term behavior of particle trajec- 
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Figure 12: Mean square displacement (r 2 (f)) vs time t for Kac diffusion 
on Sierpinski carpet (panel a) and on the deterministic cross fractal (panel 
b) at D = 1. Lines (a) to (c) refer to different values of a = 0.1 (lines 
a), a = 1 (lines b), a = 10 (lines c). Lines (d) correspond to the scaling 
( r 2 [t )) ~ f 2 . Lines (e) to the anomalous diffusion scaling ( r 2 (t )) ~ t 2 / d ™ ; 
where 2/d w = 0.94 (panel a), and 2 /d w = 0.8114 (panel b). 


tories a fractal behavior characterized by a fractal dimension different from 
3/2, that corresponds to the classical Brownian motion in Euclidean media. 

Let dx be the fractal dimension of Poisson-Kac trajectories in the long¬ 
term regime, and L(Af; t max ) the length of a finite portion of the trajectory 
in the interval t € [0, t max ) estimated using a temporal yardstick At. Hence- 
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Figure 13: Switching time probability density functions g a {r), a = 1,2, 
vs t for diffusion on the Sierpinski carpet at D = 1. Panel (a): a = 0.1, 
panel (b): a = 10. Symbols (□) refer to <?i(t), symbols (•) to <72 ('?")• The 
solid line is the bare Poisson statistics g a (r) = ae~ at , the dashed line is the 
exponential fit of the data. 


forth, for notational simplicity, we drop the reference to t max and write 
simply L(At). 

It follows from the above qualitative analysis that 

f constant At < At* 

{ At l ~ dT At > At*, 




where the crossover time At* depends on a. The fractal dimension dx is 
related to the Hurst exponent of the graph of the trajectories by the relation 

[31] 

d T = 2 - H , (41) 
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Figure 14: Switching time probability density functions <7 q ( t ), a = 1,2, vs 
r for diffusion on the deterministic cross fractal at D = 1, a = 1. Symbols 
(□) refer to gi(r), symbols (•) to <72 (t - )- The solid line is the bare Poisson 
statistics g a (r) = ae~ at , the dashed line is the exponential fit of the data. 

and, in turn, it is related to the walk dimension d w by H = l/d w , so that 

d T = 2 - • (42) 

Consequently, the length-yardstick analysis, typical of the geometric char¬ 
acterization of fractal curves EH, provides an alternative way to estimate 
d w from a single realization of the process. 

Figure [T7] depicts the graphs of L(At) vs At for several Poisson-Kac 
processes on the Siepinski carpet and on the deterministic cross fractal at 
different values of a. The realization of the process has been obtained by 
integrating eq. (fTBll with a time step h t = 2 x 10“ 3 up to t max = 4 x 10 2 , 
starting from a point inside the fractal structure chosen at random. 

The data depicted in figure [T7] agree with the scaling (14011 . For small 
At < At*, where At* decreases with cl, T(At) saturate, due to the rectifiable 
nature of the trajectories. For At > At* the fractal scaling emerges as a 
long-term property of the process. From eqs. (l40l) - ([42l) it follows that 

L(At) ~ At l/d ™- 1 , At > At* , (43) 

and the value of d w estimated from the long-term fractal scaling of L(At) 
using eq. (14311 perfectly agrees quantitatively with the value of the walk 
dimension obtained from the analysis of particle mean square displacement. 
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Figure 15: 4> co \\ vs a for Kac diffusion on fractals at D = 1. Solid lines 
represent the scaling </> co ii(a) ~ °■ Line (a) and (•) refer to the Sierpinski 
carpet (£ = 0.72), line (b) and (□) to the deterministic cross fractal (£ = 
0.79). 

As expected, the fractal scaling appears more neatly (i.e., over a larger 
range of At) if one considers higher values of a (curves (c) in figure [17)1 . 
since, for fixed D, the characteristic stochastic velocity b = y/2 D a is higher 
and the collisions with the fractal boundary, determining the occurrence of 
the emergent fractal behavior of the trajectories, stabilize the process at 
smaller At. 

6 Concluding remarks 

The statistics of Poisson-Kac processes in closed systems is profoundly af¬ 
fected by reflection conditions at the boundaries. This effect influences either 
the statistics of the switching times and the emerging long-term properties of 
the transport process. The study of Poisson-Kac diffusion on fractals clearly 
shows that the anomalies in diffusion occur as emerging properties of a local, 
almost everywhere differentiable, stochastic motion at microscales, deriving 
from the complex collision process occurring at the boundary of a fractal 
set, and determining the transition from a ballistic motion to anomalous 
Einstein scaling characterized by a walk dimension greater than 2. 

The crossover in the mean square displacement has its dynamic coun¬ 
terpart in the properties of the trajectories deriving from the analysis of 
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Figure 16: x-component x(t) of a diffusive Kac trajectory vs t on the de¬ 
terministic cross fractal at D = 1. Line (a) refers to a = 0.1, line (b) to 
a = 10. 

the scaling of the length L(At) vs the timescale At. The trajectories of 
Poisson-Kac diffusive particles admit a transition from a local smooth be¬ 
havior at shorter timescales to the emerging fractal properties characterized 
by a trajectory fractal dimension (It = 2 — l/d w . 

Two observations deserve further attention. Poisson-Kac processes ad¬ 
dressed in Section [5] are indeed useful tools to investigate, using off-lattice 
algorithms, complex transport processes in fractal and disordered media, ei¬ 
ther in the case of pure diffusive motion (as in the present work) or including 
the effects of deterministic velocity fields and potentials. At relatively small 
values of a and b, they represent the stochastic transport of a kind of a vis¬ 
coelastic phase with memory, and memory effects influence the short-term 
behavior. Conversely the long-term, long-distance properties of Poisson-Kac 
processes are those of classical Brownian motions. 

However, Poisson-Kac processes do not share with the classical Brownian 
motion all the conceptual problems associated with an infinite propagation 
velocity, and therefore they can be consistently applied also to relativistic 
problems. 

The second observation is purely technical. In order to highlight the 
long-term fractality of Poisson-Kac trajectories emerging from locally dif¬ 
ferentiable stochastic motion we have used the length/time-interval scaling 
(figure ITTjl . out of which the Holder exponent H of the stochastic trajec- 
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Figure 17: Scaling of the length L(At) vs the temporal yardstick size At on a 
diffusive Kac trajectory at D = 1. Panel (a) refers to the deterministic cross 
fractal, panel (b) to the Sierpinski carpet. Lines (a) and (□) correspond to 
a = 0.1, lines (b) and (o) to a = 1, lines (c) and (•) to a = 10. Lines (d) 
corresponds to the theoretical prediction L(At) ~ At 1 /^ -1 , where d w = 
log5/log 3 + 1 in panel (a), and d w = 2.4648 in panel (b). 


tories can be estimated. Since H is related to the walk dimension d w , the 
scaling of L(At) vs At provides a reliable, accurate and efficient way to ana- 
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lyze anomalous transport properties in fractal media, alternative to the more 
classical analysis of the scaling of the mean square displacement (?’ 2 (t)), just 
using a single particle trajectory. 
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